device,decomposed=0
loadct,5
default,field,'u'
default,extension,'xy'

eu_read_param,o=o
n = o.n
m = o.m
l = o.l
nt = o.nt
nslice = o.nslice

theta = indgen(m)*!pi/(m-1)
phi = indgen(n)*2.*!pi/(n-1)
rds = 4.88e8 ; base of the convection zone in [m/s]
rr = 6.97e8  ; Solar radii in [m/s]
r = (rds + indgen(l)*o.dz00)/rr

x = r#sin(theta)
y = r#cos(theta)

th = 180.*theta/!pi - 90.
ph = 180.*phi/!pi - 90
     
;    contour,(wa),ph,th,/fi,nl=64,/over

if (extension EQ 'xy') then begin
  filename = field+extension+'.dat'
  arr = dblarr(n,m)
endif

if (extension EQ 'yz') then begin
  filename = field+extension+'.dat'
  arr = dblarr(m,l)
endif


if (extension EQ 'xz') then begin
  filename = field+extension+'.dat'
  arr = dblarr(n,l)
endif

t = 0.D
openr,1,filename,/f77
ns = 0L
while (not eof(1)) do begin
  readu,1,arr,t
  if (ns LE 1) then begin
    aarr = arr 
    print,'Reading file: ',filename
  endif else begin
    aarr = [[[aarr]],[[arr]]]
  endelse
  wait,0.
  ns ++
end
close,1

min = min(aarr)
max = max(aarr)
print,min,max

lev = 2.*max(abs(aarr))*indgen(nlev)/float(nlev-1)-max(abs(aarr))
if (extension EQ 'xy') then begin
;  map_set,/mollweide,/grid,/noborder,/isotropic,0,80,title='!8u!6'
  for i = 0, nt/nslice-2 do begin
    print,i
    contour,smooth(aarr[*,*,i],5),ph,th,/fi,/over
    wait,0.
  endfor 
endif

if (extension EQ 'yz') then begin
  for i = 0, nt/nslice-2 do begin
    contour,smooth(transpose(aarr[*,*,i]),1),x,y,/fi,/iso,levels=lev,/over
    colorbar, pos=[.6,.1,.62,.9], div=4, range=[min(lev),max(lev)], /vert, $
          format='(F8.2)', charsize=1.5,title=field
    wait,0.1
;    contour,transpose(total(aarr[*,*,148:248],3)/100),x,y,/fi,levels=lev,/iso
    colorbar, pos=[.6,.1,.62,.9], div=4, range=[min(lev),max(lev)], /vert, $
          format='(F5.2)', charsize=1.5,title=field
  endfor 
endif

end
